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LONG  TERM  GOALS 

The  long  term  goal  of  this  research  is  to  significantly  advance  computational  methods  for  multi-scale 
flow  physics  in  geometrically  and/or  hydrodynamically  complex  oceanic  and  coastal  ocean 
environments,  through  refinements  in  the  defined  physics,  domain  definition  and  computational  grid 
resolution.  The  particular  focus  is  on  improved  coupling  of  wind  generated  short  wave  and  circulation 
models  within  the  framework  of  adaptive,  unstructured  grid  models  for  oceanic  and  coastal  waters. 

OBJECTIVES 

The  objective  of  this  project  is  the  dynamic  coupling  of  the  ADCIRC  circulation  model  and  the  short 
wind  wave  model  SWAN  using  unstructured  adaptive  computational  meshes.  In  recent  experience 
with  coupling  the  continuous  Galerkin  based  version  of  ADCIRC  with  short  wind  wave  models,  it  has 
been  found  that  wave  transformation  zones  require  high  levels  of  resolution  in  order  to  correctly 
capture  the  wave  radiation  stress  forcing  to  the  circulation  model.  However  the  wave  transformation 
zone  tends  to  shift  depending  on  the  direction  and  period  of  the  waves.  Thus  /^-adaptive  discontinuous 
Galerkin  (DG)  based  solutions  are  being  developed  to  optimize  the  application  of  high  resolution  zones 
to  correctly  capture  the  wave  coupling  within  these  zones  without  over-resolving  adjacent  areas. 
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APPROACH 


We  have  focused  on  three  areas  to  accomplish  the  defined  objectives.  The  first  area  is  implementing 
and  evaluating  /ip-adaptive  schemes  for  DG  solutions  to  the  shallow  water  equations  (SWE)  in 
anticipation  of  coupling  with  a  mesh  adaptive  version  of  the  SWAN  wave  model  which  is  under 
development.  This  will  improve  the  unstructured  grid  capabilities  for  the  ADCIRC  circulation  model 
by  implementing  Zip-adaptive  methods  to  ensure  that  sufficient  grid  resolution  is  provided  for  all 
relevant  flow  scales.  The  present  implementation  of  ADCIRC  is  based  on  both  Continuous  Galerkin 
(CG)  and  DG  finite  element  methods.  Currently  we  are  intensively  developing  the  DG  based 
algorithms.  DG  algorithms  are  particularly  well-suited  for  both  propagation  and  advection  dominated 
problems  with  or  without  sharp  gradients  in  the  forcing  function,  bathymetry,  and/or  flow.  DG 
methods  inherently  preserve  mass  perfectly  on  an  elemental  level,  which  make  them  ideal  for  coupling 
flow  and  transport  models.  Example  problems  include  flows  with  strong  eddies  such  as  those  issuing 
from  and  through  inlets;  flows  through  rapidly  varying  bathymetry  such  as  canyons  or  deep  scour 
holes  in  inlets;  and  flows  with  sharp  fronts  such  as  tidal  bores,  ebb  tide  -  waves  interactions  in  inlets 
and  wetting  fronts.  DG  methods  are  conceptually  similar  to  finite  volume  methods  although  DG 
methods  are  readily  implemented  with  higher-order  bases.  DG  methods  are  also  ideal  since  non- 
conforming  h  and  p  refinement  and  adaptivity  can  be  implemented.  In  addition,  DG  methods  are 
significantly  more  accurate  on  a  per  degree  of  freedom  basis  than  CG  methods.  Finally  DG  based 
methods  require  less  message  passing  within  a  parallel  framework  and  are  therefore  more  suitable  for 
massively  parallel  processing. 

The  second  area  is  evaluating  the  influence  of  grid  resolution  in  the  SWE  circulation  model  in  relation 
to  the  wind- wave  model  in  the  computed  wave  radiation  stress  set  up.  In  addition,  we  are  studying  the 
sequencing  of  interpolation  and  the  wave  radiation  stress  gradient  computation  and  the  influence  of 
both  the  wave  and  surge  model  grids  in  this  context. 

The  third  area  is  to  improve  the  wave  -  circulation  model  coupling  which  was  bottle  necked  by  our 
single  processor  interface  for  communication  between  the  two  models.  We  are  designing  all  paradigms 
for  very  large  scale  parallelism  on  distributed  memory  machines,  anticipating  thousands  of  processors. 
We  are  currently  using  the  globalized  output  to  localized  input/interpolation  to  communicate  between 
the  wind  wave  and  circulation  models.  We  are  working  on  a  synchronous  parallel/parallel  interface. 

WORK  COMPLETED 

Work  has  continued  on  the  DG  modeling  system  development  focusing  on  system  validation, 
improvements  in  the  automated  p-adaptivity,  implementation  and  testing  automated  ^-adaptivity, 
improvements  in  parallel  coding  and  scaling,  and  implementing  three  dimensional  baroclinic  terms. 

System  validation  work  has  focused  on  basin  scale  to  inlet  scale  problems  to  ensure  that  the  DG  based 
algorithms  are  robust  over  the  entire  range  of  momentum  balances  and  allow  the  use  of  large  scale 
domains.  Large  scale  domains  are  attractive  since  they  significantly  improve  the  coupling  mechanisms 
for  deep  ocean  water  to  propagate  onto  the  shelf  and  into  inlets  and  adjacent  floodplains  for  all 
energetic  scales  of  motion.  Validation  exercises  have  focused  on  tides.  Tide  simulations  include  tidal 
boundary  forcing  functions  and  tidal  potential  functions. 

DG  development  work  has  advanced  significantly.  Improvements  have  been  made  in  the  robustness 
and  accuracy  of  the  DG  based  wetting/drying  algorithm,  including  proving  convergence  for  a  range  of 
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problems.  The  DG  code  now  features  both  p-  (polynomial  order)  and  standard  h-  (grid) 
refinement/adaptivity  options.  The  /^-adaptivity  is  made  simple  with  the  use  of  a  hierarchical  basis  and 
refinement  is  simply  based  on  adding  higher-order  basis  functions.  Derefinement  is  implemented  by 
discarding  higher-order  basis  functions,  equivalent  to  an  L2  projection  of  solution  onto  lower  order 
space.  A  "sensor"  proposed  by  Burbeau  and  Sagaut  (2005)  is  used  to  decide  where  to  p- refine.  The  h 
data  structure  of  Demkowicz  (2006)  has  been  implemented.  This  structure  stores  only  minimal  mesh 
information  and  no  information  for  new  elements  is  stored  explicitly  whatsoever.  Genealogical 
information  stored  for  element  nodes  -  i.e.  children  and  parents  of  nodes  and  element  connectivities 
are  then  recovered  from  this  information  using  fast  data  structure  supporting  algorithms.  The  h- 
adaptive  algorithm  performs  standard  h4  isometric  refinement  (Figure  1).  The  method  allows  the  use 
of  hanging  nodes  but  is  limited  to  so-called  1-irregular  meshes  (Figure  2).  A  discontinuity  detector  is 
used  to  decide  where  to  h- refine  (Krivodonova  et  al.  2004).  Note  both  the  indicators  of  h-  and  p- 
adaptation  are  based  on  quantities  that  are  easy  to  compute  or  already  available  and  maintain  the  local 
character  of  the  DG  method. 

Improvements  have  been  made  in  both  the  efficiency  of  the  DG  code  as  well  as  the  parallel 
implementation.  Through  loop  level  optimization,  we  have  improved  the  efficiency  of  the  DG  code  so 
that  on  a  per  degree  of  freedom  basis  it  is  only  3.5  times  more  expensive  to  run  than  the  CG  code 
(compared  to  8  times  more  expensive  prior  to  optimization).  It  is  noted  that  the  DG  code  is 
significantly  more  accurate  than  CG  on  a  per  degree  of  freedom  basis.  Parallelization  has  also  been 
improved  so  that  linear  speedups  are  maintained  much  longer  than  CG  based  methods.  Parallelization 
is  achieved  through  domain  decomposition  with  a  shared  overlapping  layer  of  element  with  inter¬ 
processor  communications  based  on  embedded  shared  layer  elements  on  one  processor  sending 
polynomial  expansion  coefficients  to  the  corresponding  exterior  element  in  the  adjacent  processor. 

This  allows  all  computations  including  the  Roe  flux  computation  to  be  performed  on  each  individual 
processor. 

Since  the  goal  is  to  develop  multi-scale  adaptive  flow  computation  capabilities  for  the  entire  range  of 
coastal  ocean  physics  and  problems,  we  have  also  investigated  the  application  of  DG  methods  to  three- 
dimensional,  baroclinic  models.  A  three-dimensional  DG  model  was  developed  by  Dawson  and 
Aizinger  (2007),  and  has  been  extended  from  barotropic  to  baroclinic  mode.  The  model  uses  a  z-grid 
with  general  prismatic  elements  and  has  the  capability  to  follow  a  sloped  (piecewise  linear)  bottom 
bed.  The  method  is  based  on  the  use  of  discontinuous,  piecewise  polynomial  approximating  functions 
for  each  primary  variable,  defined  over  each  element.  The  primary  focus  during  the  past  year  was  the 
incorporation  of  various  turbulence  closure  models  into  the  code  and  the  preliminary  testing  of  these 
models.  Turbulence  models  that  were  implemented  include  two  algebraic  (zeroth  order)  models,  a  one 
equation  model,  and  various  two  equation  models.  Results  using  different  orders  of  approximating 
spaces,  from  piecewise  constant  to  piecewise  quadratic,  indicate  that  piecewise  quadratic 
approximations  of  kinetic  energy  are  superior  at  resolving  the  boundary  layer  at  the  sea  bed. 

We  have  investigated  vertical  eddy  viscosity  models  of  various  levels  of  computational  and  conceptual 
complexity.  In  order  of  increasing  complexity  those  include  a  constant  eddy  viscosity  coefficient,  an 
algebraic  (zeroth  order)  model,  as  well  as  one  and  two  equation  models.  The  simplest  model  amounts 
to  explicitly  specifying  diagonal  entries  to  the  tensors  of  eddy  viscosity/diffusivity  coefficients  for  all 
variables.  Two  algebraic  models  implemented  in  our  DG  code  are  due  to  Davies  (1986)  and  give  good 
results  at  a  reasonable  computational  cost  in  cases  where  accurate  vertical  resolution  of  flow  is  not 
important.  The  first  order  vertical  eddy  viscosity  closure  model  solves  a  transport  equation  for  the 
turbulent  kinetic  energy  in  addition  to  the  mass,  momentum,  and  species  transport  equations.  The 
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turbulent  mixing  length  is  computed  algebraically  in  this  model  (see  Delft3D-Flow  manual,  2005).  The 
second  order  closure  model  we  implemented  is  based  on  the  generic  turbulence  length  scale  model 
proposed  by  Warner  et  al.  (2005).  The  main  advantage  of  this  formulation  is  the  ability  to  switch 
between  several  two  equation  models,  including  k-epsilon  and  Mellor-Yamada,  by  changing  a  few 
constant  parameters. 

In  addition  to  developing  robust  basin  to  channel  scale  DG  based  /^-adaptive  algorithms,  we  are 
studying  the  effects  of  wave  radiation  stress  induced  set  up  in  hurricane  environments.  Using  models 
developed  for  the  USACE  in  Southern  Louisiana,  we  are  examining  the  degree  of  wave  set  up.  We  are 
designing  ^-adaptivity  in  order  to  minimize  the  high  level  of  resolution  in  the  highly  refined  broad 
wave  transformation  zones  that  we  have  built  into  the  current  static  models. 

In  the  area  of  inter-model  communications  we  have  refined  the  global  to  local  paradigms  and 
interpolation.  We  have  also  improved  code  efficiency  by  specializing  output  processors. 

RESULTS 

The  DG  based  code  has  been  applied  to  a  domain  encompassing  the  western  north  Atlantic  and  the 
entire  Gulf  of  Mexico  and  Caribbean  Sea.  This  domain  allows  tides  to  interact  between  basins,  be 
amplified  within  the  resonant  Gulf  and  Caribbean  and  propagate  onto  the  shelves.  The  coarse  grid 
applied,  shown  in  Figure  3,  works  very  well  with  the  DG  solution  as  is  illustrated  in  the  sample  tidal 
comparisons  shown  in  Figures  3  and  4.  A  much  more  nonlinear  case  is  the  tide  computation  in  Willapa 
Bay  shown  in  Figure  5.  This  case  tests  the  advective  terms  as  well  as  the  new  DG  wetting/drying 
algorithm  since  the  tide  range  is  large  (more  than  4  m),  there  is  a  narrow  inlet  to  a  large  back  bay  and 
there  are  large  tidal  flats  with  a  significant  portion  of  the  bay  drying  during  the  tidal  cycle.  This 
problem  encompasses  basin  to  channel  scales.  Results  at  NOS  stations  are  compared  to  data  in  Figure 
6. 

The  /^-adaptive  solutions  have  now  been  applied  to  Grays  Harbor,  a  highly  advection  dominated  case 
characterized  by  a  deep  inlet  channel  and  long  jetties  which  are  partially  submerged.  The  grid  and 
bathymetric  detail  is  shown  in  Figures  7  and  8.  /7-adaptivity  is  shown  in  Figures  9  and  1 1  for  flood  and 
ebb  tides  respectively,  while  the  corresponding  currents  are  shown  in  Figures  10  and  12.  We  note  that 
eddies  as  well  as  the  strong  shear  are  nicely  picked  up  by  the  p  refinement  sensors. 

The  /7-adaptive  algorithms  have  been  tested  on  the  idealized  dam  break  problem  illustrated  in  Figure 
13.  The  colors  indicate  surface  water  elevation.  Initially  the  left  and  right  sides  of  the  basin  are  at 
different  water  levels  and  a  sluice  is  then  instantly  opened  between  the  two  portions  of  the  basin.  It  is 
noted  that  the  shock  propagates  from  the  left  side  to  the  right  side  and  as  it  does  so  the  /7-adaptivity 
refines  the  grid  to  capture  the  large  gradients. 

To  evaluate  the  parallel  performance  of  the  DG  model,  we  tested  a  medium  size  problem  and 
compared  the  performance  to  that  of  our  CG  based  code.  The  problem  consists  of  the  western  North 
Atlantic  Ocean,  Gulf  of  Mexico  and  Caribbean  Sea  and  in  the  grid  contains  492,179  elements  and 
254,565  nodes.  Simulations  were  performed  on  the  TACC’s  Lonestar  Dell  Linux  Cluster,  which 
consists  of  5,840  cores  within  1,460  Dell  PowerEdge  1955  blades  each  with  two  Xeon  5100  series 
2.66GHz  dual-core  processors.  A  simulation  was  performed  on  16,  32,  64,  128,  256,  512,  and  1024 
cores.  Figure  14  shows  the  speedup,  defined  as  the  ratio  of  the  wall  clock  time  of  a  single  core  run  over 
the  wall  clock  time  of  a  parallel  run.  For  this  test  case,  the  ADCIRC-CG  code  exhibits  linear  speed-up 
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to  64  cores  but  tails  off  beyond  128  cores.  This  behavior  is  expected  for  problems  of  fixed  size;  our 
experience  to  date  indicates  that  the  code  begins  to  lose  parallel  efficiency  at  around  2000  finite 
element  mesh  nodes/core.  The  ADCIRC-DG  code  scales  close  to  linearly  up  to  256  processors  and 
continues  to  scale  with  efficiency  above  50%  to  1024  cores.  The  corresponding  efficiency  of  the  CG 
and  DG  codes  is  shown  in  Figure  15.  It  is  noted  that  more  parallel  efficiency  improvements  are 
possible  for  the  DG  code  and  that  scaling  and  efficiency  improve  for  larger  problems. 

As  a  test  run  for  the  three  dimensional  implementation,  we  simulate  tide-driven  flow  in  the  Bight  of 
Abaco.  The  domain  bathymetry  and  the  coarse  mesh  are  shown  in  Figures  16  and  17.  A  standard  tidal 
forcing  of  five  constituents  was  imposed  at  the  open  sea  boundary.  In  Figures  18  through  22,  we 
compare  the  free  surface  elevation  computed  using  piecewise  linear  approximation  spaces  and  up  to  5 
layers  in  the  vertical  direction,  using  algebraic  (Figure  18),  Delft3D  order  1  (Figure  19),  k-epsilon 
(Figure  20),  k-omega  (Figure  21)  and  the  "generic"  order  2  model  (Figure  22).  The  results,  including 
those  obtained  using  the  algebraic  turbulence  closure  model,  are  closely  matched. 

For  studying  wave  radiation  stress  gradient  induced  set  up,  we  have  been  applying  the  SL15 
domain/grid  for  Southern  Louisiana  shown  in  Figures  23-27.  This  grid  was  developed  for  the 
USACE/FEMA  Joint  Coastal  Surge  Study  and  has  very  high  levels  of  shelf,  coastal  and  inland 
resolution  in  order  to  compute  hurricane  storm  surge.  In  addition,  as  can  be  seen  in  Figure  25,  a  wide 
swath  of  high  resolution  (80m  and  smaller)  has  been  specified  in  the  coastal  wave  transformation  zone. 
In  fact,  about  40%  of  the  2.17  million  computational  points  resolve  the  wave  transformation  zone.  The 
corresponding  wave  radiation  stress  high  gradient  zones  for  given  storms  and  parts  of  a  specific  storm 
are  much  narrower  than  the  swaths  of  high  resolution,  although  the  zones  shift  in  time  with  the  storm 
and  wave  period.  We  have  been  working  on  the  implementation  of  the  DG  code  and  the  /^-adaptive 
DG  code  to  this  problem.  A  CG  evaluation  of  wave  induced  radiation  stress  for  Hurricanes  Katrina  and 
Rita  indicates  that  wave  radiation  stress  induced  set  up  along  the  coast  is  generally  less  than  0.5  m  and 
up  to  1  m  in  well  focused  areas  that  do  not  have  extensive  back-bay  areas  for  the  water  to  spread  out. 
Wave  radiation  induced  set  up  is  less  for  Rita  than  for  Katrina  because  of  the  predominant  winds 
coming  off  of  land  during  much  of  the  storm.  We  are  systematically  going  through  grid  convergence 
studies  to  define  levels  of  required  resolution. 

IMPACT/APPLICATIONS 

This  work  will  significantly  improve  the  accuracy  of  the  computed  physics  for  coastal  ocean 
hydrodynamic  computations,  especially  for  multi-scale  physics  since  the  DG  based  adaptivity  will 
automatically  and  optimally  sense  where  the  mesh  refinement  is  needed.  It  will  also  reduce  model 
development  times  since  base  meshes  only  need  to  represent  the  geometry  and  not  the  hydrodynamics. 
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Figure  2:  1-irregular  meshes  with  hanging  nodes . 
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Figure  3:  Coarse  grid  Western  North  Atlantic  domain  with  sample  tide 
comparisons  to  Fourier  re-synthesized  measured  tidal  data. 
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Figure  4:  Sample  tide  comparisons  to  Fourier  re-synthesized  measured 
tidal  data  along  U.S.  coastline. 
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Figure  5:  Willapa  Bay  grid  and  topography/detail  within  the  bay 
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Willapa  Bay,  WA 


Figure  6:  Willapa  Bay  tidal  validation  comparing  model  tides  and 
Fourier  re-synthesized  measured  tidal  data. 
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Figure  7:  Grays  Harbor  p-adaptivity  testing  to  resolve  eddy  features  around  jetties. 


Figure  8:  Grays  Harbor  grid  and  bathymetric  detail. 
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Figure  9:  Grays  Harbor  p-adaptivity  showing  p=3  elements  in  purple 
and  p=l  elements  in  blue  on  the  flood  tide. 


Figure  1 0:  Grays  Harbor  p-adaptivity  showing  surface  water  elevations 
and  current  vectors  on  the  flood  tide. 
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Figure  11:  Grays  Harbor  p-adaptivity  showing  p=3  elements  in  purple 
and  p=l  elements  in  blue  on  the  ebb  tide. 


Figure  12:  Grays  Harbor  p-adaptivity  showing  surface  water  elevations 
and  current  vectors  on  the  ebb  tide 
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Figure  13:  Idealized  dam  break  problem  with  levels  of  h-adaptivity 
and  colors  indicating  water  level. 


Speedup  =  (Single  processor  wall  clock  time)/(Parallel  wall  clock  time) 


Figure  14:  Idealized  parallel  speedup  of  the  ADCIRC-CG  and  ADCIRC-DG 
codes  up  to  1024  cores  for  a  492,1 79  element  model. 
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Efficiency 


Efficiency  =  (Speedup)/(#  of  processors) 


Figure  15:  Parallel  efficiency  of  the  ADCIRC-CG  and  ADCIRC-DG  codes 
up  to  1024  cores  for  a  492,1 79  element  model. 
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Figure  16:  Bathymetry  of  the  Bight  of  Abaco,  Bahamas. 
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P.w.  li nears,  5  layers,  algebraic  closure  scheme 
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Figure  18:  Free  surface  elevation  at  t=l,000,000  seconds  in  the 
Bight  of  Abaco  using  algebraic  closure. 


jP.w.  li nears,  5  layers,  DeJft3D  scheme! 


Figure  19:  Free  surface  elevation  at  t=l,000,000  seconds  in  the 
Bight  of  Abaco  using  Delft  3D  order  1  closure. 


16 


P  .w.  linears,  5  layers,  k-ep  s  i  I  on  sche  me 
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Figure  20:  Free  surface  elevation  at  t=l,000,000  seconds  in  the 
Bight  of  Abaco  using  k-epsilon  closure. 


Figure  21:  Free  surface  elevation  at  t=l,000,000  seconds  in  the 
Bight  of  Abaco  using  k-omega  closure. 
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P.w.  linears^S  layers,  generic  scheme 
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Figure  22:  Free  surface  elevation  at  t=l,000,000  seconds  in  the 
Bight  of  Abaco  using  the  generic  order  2  closure . 
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Figure  23:  SL15  Computational  domain  and  finite  element  grid  for  studying  riverine,  tidally,  storm 
surge  and  wave  radiation  stress  gradient  induced  flow  in  Louisiana  and  Mississippi. 
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Figure  24:  SL15  bathymetry  for  studying  studying  riverine,  tidally,  storm  surge  and 
wave  radiation  stress  gradient  induced  flow  in  Louisiana  and  Mississippi. 


Figure  25:  SL15  Computational  grid  detail  for  studying  riverine,  tidally,  storm  surge  and  wave 
radiation  stress  gradient  induced  flow  in  eastern  Louisiana  and  western  Mississippi 
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Figure  26:  SL15  bathymetry  and  topography  for  studying  riverine,  tidally,  storm  surge  and  wave 
radiation  stress  gradient  induced  flow  in  eastern  Louisiana  and  western  Mississippi 


Figure  27:  Computed  wave  radiation  stress  induced  set  up  for  eastern 
Louisiana  and  western  Mississippi  during  Hurricane  Katrina. 
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Figure  28:  Computed  wave  radiation  stress  induced  set  up  for  eastern  Louisiana  and  western 

Mississippi  during  Hurricane  Rita. 


Figure  29:  Computed  wave  radiation  stress  induced  set  up  for  western 
Louisiana  during  Hurricane  Rita. 
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